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Abstract 



Early work extending the Kohn-Sham theory to excited states was based 
on replacing the study of the ground-state energy as a functional of the 
ground-state density by a study of an ensemble average of the Hamiltonian as 
a functional of the corresponding average density. We suggest and develop an 
alternative to this description of excited states that utilizes the matrix of the 
density operator taken between any two states of the included space. Such 
an approach provides more detailed information about the states included, 
for example, transition probabilities between discrete states of local one-body 
operators. The new theory is also based on a variational principle for the 
trace of the Hamiltonian over the space of states that we wish to describe 
viewed, however, as a functional of the associated array of matrix elements 
of the density. It finds expression in a matrix version of Kohn-Sham theory. 
To illustrate the formalism, we study a suitably defined weak-coupling limit, 
which is our equivalent of the linear response approximation. On this ba- 
sis, we derive an eigenvalue equation that has the same form as an equation 
derived directly from the time-dependent Kohn-Sham equation and applied 
recently with considerable success to molecular excitations. We provide an 
independent proof, within the defined approximations, that the eigenvalues 
can be interpreted as true excitation energies, 
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I. INTRODUCTION 



Density functional theory (DFT) was designed originally as a theory of the ground-state 
of a many-particle system For an extension to include the calculation of excitation 

energies, several lines of thought have been developed. The earliest was based on a minimum 
principle ||[7j for the trace of the Hamiltonian over a set of the lowest-energy eigenstates 
of the system. This theory has been extended to a suitably weighted sum over the same 
set of eigenstates M. The expanded version of the Hohenberg-Kohn theorem, in either 
case, is that the average energy is a unique functional of the corresponding average density. 
Excitation energies are (essentially) obtained by taking differences between averages over 
almost overlapping sets. This approach has not been developed beyond the cited work. 

Recently, considerable attention has been focused on the development of other methods 
for studying excitation energies. One approach is based on time-dependent density functional 
theory (TDDFT) Here one studies the linear response of the time-dependent density 

to a time-dependent external field. The Fourier transform of the susceptibility (density- 
density correlation function), which is the essential ingredient for the calculation of dynamic 
polarizabilities, has poles at the true eigenstates of the system. 

From these elements two different formalisms for the calculation of energies and po- 
larizabilities have been derived, termed density based and density-matrix based coupled 
Kohn-Sham schemes (CKS) in a recent publication fT5| . Both methods have been applied 



successfully, the density method in (|Tt],|nj, the density matrix method more extensively in 
[|T4^f2~T||. The mathematical equivalence of the two formalisms derived from TDDFT has 
been established in Ref . |15| , so that it is only a matter of numerical convenience which one 



is used in practice. To argue the utility of the formalism developed in this paper, we shall 
show that in a suitable approximation it yields the eigenvalue equation of the density-matrix 
CKS scheme. It also permits the calculation of electromagnetic transition rates between the 
excited states studied and the ground state without the explicit use of wave functions. 

We also call attention to several recent studies of the excited state problem that involve 
extensions of the variationally based KS theory to individual excited states [p2,p3|. For 



these methods, as well, successful applications have been made to simple systems. Improved 
exchange and correlation kernels necessary for all these methods and a connection with 
many-body perturbation theory are discussed in [ 23 1 , whereas in [21j] an exchange-correlation 



potential is suggested to provide more accurate continuum KS orbitals needed for excited 
state and polarizability calculations. 

In this paper, we return to a study of the trace variational principle p6|-p8|j . Instead 
of considering the average energy as a functional of the average density, however, we argue 
for the introduction of a matrix array of densities, i. e., all matrix elements of the density 
operator among all states of a chosen ensemble, and for an investigation of the average 
energy as a functional of this matrix array. In Sec. II we present arguments to indicate 
how the Hohenberg-Kohn (HK) analysis can be extended to this case yielding a variational 
equation for the matrix. We subsequently (Sec. Ill) generalize the KS analysis, deriving 
a matrix Kohn-Sham equation (MKS), that contains not only the expected ingredient, a 
matrix effective potential, but also a matrix of Lagrange multipliers arising from number 
conservation in each state of the chosen subset; this matrix can be diagonalized, but not 
otherwise transformed away. By combining solutions of the MKS equations, we can construct 
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the density array. 

As an application of this theory, we study, in Sec. IV, the MKS equations in the weak- 
coupling limit, which is the equivalent, in our approach, to linear response theory. In this 
limit, we include only the ground state and excited states characterized (largely) by one 
'quasiparticle-quasihole' excitations of the ground-state. Higher excited states are incorpo- 
rated via simple assumptions concerning their properties. The major result of this analysis 
is an eigenvalue equation for the aforementioned Lagrange multipliers (relative to their 
ground-state value) that has the same form as the eigenvalue equation of the density-matrix 
based CKS. Assuming that the ground-state KS problem has been solved, the major un- 
known ingredient in these equations, an exchange-correlation interaction, can be identified 
with the corresponding quantity utilized in current applications, which in practice require 
an adiabatic approximation. 

Comparison of our result with the existing formalism establishes the identity of the 
eigenvalues of our equation with excitation energies of the system. These energies can be 
calculated as well from a difference of adjacent averages of the Hamiltonian. By this means, 
in Sec. V, we establish within the framework of our formalism the physical meaning of the 
eigenvalues. In a concluding section, we summarize our considerations. 



II. HOHENBERG-KOHN ARGUMENTS 

The Hamiltonian is written as 

H = f + V + W + Y, (2.1) 

the sum of the kinetic energy, the electrostatic interaction of the electrons with the nucleus, 
the Coulomb repulsion of the electrons, and an additional fictitious external source term 
that will be set to zero for actual calculations. The various terms have the forms (x stands 
for the space-spin pair (r, s)), in atomic units, 

dx^( x )(_Iv 2 Mx) 

$t4, (2.2) 

V = J rfx^ f (x)^(x)t;(r), (2.3) 

W = [ dxdx! -. 1 > t (x)^ t (x / )^(x / )^(x), (2.4) 
J |r — r'| 

Y = J dxrfx'y(x,x')r}(x,x'), (2.5) 

fj = ^ t (x)^(x)^ t (x')^(x / ). (2.6) 

The interaction term Y is a combination of one and two body forces. It serves as a device to 
establish the dependence of the theory on the extended density defined in Eq. ( |2.18| ), but is 
set equal to zero thereafter. This procedure is akin to the use of an external magnetic field 
in standard density functional theory in order to exhibit the possible dependence on spin 
polarization. For the traces of these operators over the ensembles introduced below, we use 
the same symbols without hats. 
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In the following we shall base our arguments on the variational principle for the trace 



of the Hamiltonian over the lowest M eigenstates of a many body system p-p|,^6|-p8| . We 
consider the case where the (M + l) st state has a higher energy than the Mth state, although 
this criterion is not absolutely necessary. 
Let 

S = {\I)} (2.7) 

be the space of states included (J = 1...M). For any operator O, we define the restricted 
trace 



M 

(M) = J2(I\d\I), (2.8) 

7=1 



where it is convenient in the further development not to divide by M. 

We then consider a set of propositions formulated in imitation of the Hohenberg-Kohn 
(HK) theorem 0: 



'i) Every choice of a function y(x, x') in (|2.5|) determines an M-dimensional space S 
through the solution of the Schrodinger equation. Two different functions y ^ y' will yield 
different S ^ S', provided Y [yj and Y[y'] differ by more than a diagonal matrix in the space 
S. To see this, suppose that Y[y] and Y[y'] yield the same space S. From the Schrodinger 
equations for the state \I), we obtain by subtraction that 

{Y[y]-Y[y']}\I} = (Ei-E' I )\I}, (2.9) 

where the E's are the corresponding eigenvalues. Thus, with I' ^ I, 

(I'\Y[ y ]\I) = (I'\Y[ y ']\I), (2.10) 

and 

M 

Y[y]-Y[y'] = j:(Ei-E' I )\I)(I\. (2.11) 
i=i 

Potentials that satisfy this relation will be considered equivalent. 

(ii) <S determines the correlation function ^(x, x') = ^(^^(X) x ')l-0- This relationship 
is single-valued and invertible. This can be proved by an adaptation of the standard HK 
argument, as we now show. Suppose that 

S^T], S'^S^rf. (2.12) 

It follows that 7] 7^ 7]' . We prove this by using the trace variational principle (valid by 
construction of the sets, as indicated above) to establish two inequalities, 

Hs[ y }<H s >[y'} + f(y-y'W, (2.13) 
H s ,[ y ']<H s [ y ] + J(y f -y)v- (2.14) 



4 



Here, for example, Hg[y] is the ensemble average of H over the set S, where it is further 
emphasized that this average is a functional of y. Adding ( |2.13 ) and ( |2.14j ) and assuming 
that r\ = rf, we obtain the usual contradiction 

H s [y] + H s ,[y'] < H s >[y'] + H s [y}. (2.15) 

Thus S is a single- valued functional of r\. 

Considering if to be a functional of rj, we write the variational principle in the form 



/x tt 
—5r] = 0. (2.16) 
01] 



We shall not attempt, however, to implement the variational principle in this version. In- 
stead, using completeness, we introduce the relation 

M oo 

r ? (x,x') = EE( / l^ t ( x M x )l // ) 

7=1 1'=l 

x(/^VMx')|/>- (2.17) 

As long as M is finite, this sum is asymmetric in the two sets of indices. However our 
aim is to utilize as variational parameters the quantities 

n(x) /a = (J|^t( x )^( x )|/') {I, I'eS}, (2.18) 

which constitute the elements of an M x M square matrix n. There are two arguments 
that can be put forward to satisfy ourselves that we can consider the correlation function 
r/(x, x') to be a functional of this more restricted set of matrix elements. The first and 
more satisfactory one at this stage is to let the number M of states included in the trace 
become large enough so that either pointwise or in the norm the fractional contribution from 
pieces connecting I < M to / > M is negligible. The second is to ask for the temporary 
indulgence of the reader, by guaranteeing to show that in any application, these matrix 
elements can either be ignored or expressed explicitly in terms of those of n. We call this 
a closure approximation. When we turn to the application developed in Sees. IV and V, 
we shall see, that the accuracy of results obtained there depend on a choice of a closure 
approximation that is historically tied to linear response theory and the associated random 
phase approximation, namely the approximate characterization of the simplest excited states 
as boson excitations. 

In consequence of the arguments given above, we replace the variational principle (|2.16|) 
by the form 

6H = ( ^-5n. (2.19) 
J on 



From Eq. (|2~T9p we can derive a variational equation by imposing the number conservation 



constraints. If N is the number of electrons, we have 

G?xn(x) 7 j, = N6 ir . (2.20) 
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Introducing a set of Lagrange multipliers pjp, we now write 

6H - mp J 6n{x)n = 0, (2.21) 



and conclude that 

SH 



Sn(x) ri 



MP ■ (2.22) 



In concluding this section, we emphasize once more that our purpose in introducing the 
unusual "external" field Y was to provide a pathway leading to the formulation, Eqs. (2.19)- 



( (2.22j ), which form the basis for the further development. These relations also apply to 
systems with Y = 0, which are the only systems that interest us henceforth. 

III. GENERALIZED KOHN-SHAM SCHEME 

n(x)/// is the limit x — > x' of the off-diagonal one-body density matrix 

p(x/|x70 = (/V t (x / )^(x)|/). (3.1) 

Since p is Hermitian and (as we shall show below) positive semi-definite in the direct product 
space labeled by (x, I), it can be brought to diagonal form, a move that generalizes the 
concept of natural orbitals. We may thus write 

p(x/|x'/') = E^(x/) < f}(xT), (3.2) 
j 

Aj > 0, (3.3) 

M , 

J2 J dx$}(x/)$^(x/) = 6 j j,, (3.4) 



i=i 



J dxp(xI\xI') = N5 IP . (3.5) 

Here Eqs. ( p.2|) and ( |3.3|) define the eigenfunctions and eigenvalues of the generalized density 
matrix, ( |3.4|) expresses the property that the $j(x/) are unit eigenvectors in the space 
labeled jointly by the single-particle coordinates and the eigenvalues of the states in the set 
iS, and ( p.5|) expresses number conservation. It follows from these equations that 

M . 

£ / dx p(x/|x/) =J2\j = NM. (3.6) 
i=i J j 

Before continuing the development, we interject the proof that p(x/|x'J / ) is positive 
semi-definite. Toward this end we compare Eq. (3^) with the form that follows directly 
from its definition 

p(x/|xT) = J2(I'\^(x')\K)(Kmx)\I}, (3.7) 

K 

where K is a complete set of intermediate states. Invoking the orthonormality relations 
flOD , we thus conclude that for eigenvalues Aj ^ 0, 
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M . 

^ = ElE/^$}(x,/)(K|^(x)|/)| 2 >0. (3.8) 

K 1=1 J 

In imitation of ground-state KS theory, we introduce a mapping from the off-diagonal 
density to a quasi-independent-particle off-diagonal density 

72(x) /7 , ^n s (x) /r , (3.9) 

n s (x) n , = ^(x/^x/'), (3.10) 

J 

J2 / dx<p*j{xl)<pj,(xl) = 5jj,, (3.11) 

J dxn s (x) ;;i = N5 ir - (3.12) 

Though we use the same symbol J to label orbitals as for the case of natural orbitals, here 
the similarity stops. For the latter, J is, in principle, an unbounded set. For the present 
alternative, the set labeled by J is strictly a finite set as determined by the sum (cf. (|3.6|) ), 

J2l=NM. (3.13) 
j 

Though the labels / have the same meaning in both the original space and in the mapped 
space, the fact that the labels J do not, leaves a vagueness at present about the significance 
of this quantum number. We shall see, however, that once we turn to applications in Sec. IV, 
there will be no problem in providing a precise characterization of this quantum label. 

We next show how the variational principle may be used to obtain equations for the 
orbitals (fj so that in fact the matrices n and n s are equal. We utilize the variational 
principle in the form 

?/ dx *Jb)^ (x/)+ c ' c -=°- (3 - 14) 



ji 



Imitating the procedure for the ground-state theory, we decompose the trace of the Hamil- 
tonian, 

H = T S + (V + W + T -T s ), (3.15) 
T S = J2[ <Pj*PJ- ( 3 - 16 ) 

Enforcing the equality of n and n s , we define an effective single-particle potential matrix, 

v s ^) ir = - /. (V + W + T- T s ), (3.17) 

(V + W + T -T s ). (3.18) 



6n s (-x) ri 



The discussion of the decomposition of this matrix single-particle operator into constituent 
parts of interest will be taken up in Sec. IV. 
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With the help of Eqs. (|3.15| - |3~T8l) , we derive from the variational principle ( |3.14j ) the 
conditions 



J2 f dx8p*j(xI)[T5 ir +w s (x) // ^ J (x/ / ) + c.c. = 0. 
j 1 1' J 



(3.19) 



To derive generalized single-particle equations of motion from the variational principle, we 
add the constraint conditions 



~ E / dxS(pj{3d)[€j6 n > + z/(x) /7 /]v?j(x/') + c.c. = 
jir J 



(3.20) 



Here ej is the Lagrange multiplier for the normalization condition contained as part of ( |3.11|) . 
(As usual, the orthogonality condition need not be imposed, since it will be automatically 
satisfied by the solutions of the emerging equations.) The unfamiliar term containing the 
Lagrange multiplier matrix ^(x)/// has the form of an additional potential matrix, whose 
purpose is to enforce the condition |HJ that n = n s . We shall study this quantity further 
below. Combining Eqs. ( |3.19|) and ( |3.20| ), we derive (together with it complex conjugate) 
the generalized single-particle equation 



ej^j(xf) = Y\ t5 h> + v s (x) n > - z/ s (x)//,]v?j(x/') 



(3.21) 



At this juncture it is appropriate to wonder how ( |3.21|) is related to the density and 
density-matrix forms of CKS theory. We cannot expect a general connection, since the latter 
describes the consequences of the application of a time-dependent external field, whereas in 
the theory under development, the "time dependence" is a purely internal matter expressed 
by an off-diagonal array of densities and effective potentials. Nevertheless, a connection 
between the two formalisms will be made for the application studied in Sec. IV, the so- 
called weak-coupling limit. In effect, it will be shown that this limit contains the same 
physics as the combination of TDDFT and linear response theory, and thus solves, in part, 
the task of establishing the "usefulness" of our approach. 

We conclude the present section by showing that (cf. Eq. (|2.22|) ) 



V-ii'i 



(3.22) 



up to a global additive constant. It is thus a non-trivial matrix and cannot be absorbed 
into the eigenvalues ej. To prove ( |3.22 ), we can work backwards from the sum of ( |3.19j ) and 
( |3.20| ) to the equation 





E 

E 

a' 

E 



<ix 



dx 



dx 



6H 



5n s (x.) ir 

SH 
5n(ic)n> 

5H 



vWir 



6n s 



i/(x) 



ii' 



5n(x.] 



i'i 



w 



5ra(x) 



pi- 



(3.23) 
(3.24) 
(3.25) 



In passing from ( |3.23p to ( |3.24| ), we have used the equality n s = n. In writing ( ^.25| ), we 
have repeated ( |2.21| ). Comparing ( |3.24| ) with ( |3.25| ), we arrive at ( |3.22|) , again up to an 
additive constant. 
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IV. APPLICATION TO THE WEAK COUPLING LIMIT 



In the course of this section, we shall transform and approximate Eq. ( p.21| ), leading to 



an eigenvalue equation that will determine off-diagonal elements of the matrix n. We shall 
do so in an approximation, the weak-coupling approximation, that is equivalent to a linear 
response approach. Assuming that the matrix fi can be chosen diagonal (see immediately 
below), the eigenvalues are the quantities 

^i = Hn-Hm- (4-1) 

The proof that the matrix fi can be chosen diagonal goes as follows: Though we trace over 
a set of states labeled / and originally identified as eigenstates of the reference system, the 
entire formalism is invariant under a unitary transformation within the included space. Such 
a transformation can be chosen to diagonalize /i if it is not already diagonal. The relation 
of the quantities in Eq. ( J4.1[ ) to the excitation energies of the system is not immediately 
apparent. The main result of this section suggests that they are equal. The proof that they 
are is given in Sec. V. 

Though the derivation of the main result of this section, the eigenvalue equation, can be 
carried out directly from the generalized KS equation, we present the discussion in a form 
that makes more immediate contact with the density functional form of the theory. The 
first step, which is completely general, is to transform Eq. (|3.21|) into an equation for the 
matrix nf 7 ,(x, x'). First rewrite Eq. ( |3.21p , remembering Eq. ( |4.1| ), as 

ej<pj(xl) = E(^/'( x ) " ujSu^jixI'), (4.2) 
v 

Recalling the definition 

r^,(xx') = £^(x/V}(xT), (4.3) 
J 

we can form from Eq. ( |4.2|) and its complex conjugate two equivalent but distinct values of 
the sum J2j e JV 9 j( x -^) ( ^j( x/ ^ / )- The difference of these forms yields the generalized density- 
matrix equation (using summation convention) 

n; J( (xx')(w/- - loj) = ^[nf J „(xx')/iJ«//(x / ) - /i* / ,,(x)nj,, J ,(xx / )], (4.4) 

/" 

that will provide the starting point for our further considerations. 
We note that by introducing time- dependent matrix elements 

Ojp(t) = Ojp exp[— i(u>i — u)ji)t], (4.5) 

where O takes on the values n s and h s , Eq. (|4.4[), may be written in the form 



-z|n s (t) = [n s (t),h s (t)]. (4.6) 

This resembles the fundamental equation of TDDFT, in density matrix form, except that 
the bold-face type reminds us that we are dealing with quantum-mechanical operators rather 
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than c-numbers. This can be converted into a form of TDDFT, however, by assuming the 
existence of a wave packet that is a linear combination of the ground state and excited 
states of interest, for which we can also replace the average of the products that appear in 
the commutator by the product of the averages. 

In the following we concentrate on the study of the weak coupling approximation to 
Eq. (|4.4j) . In this approximation one imagines that the set of states ( |2.7|) can be divided into 
the ground-state 1 = and states characterized as dominant v quasiparticle- v quasihole 
excitations (denoted by I u ) with respect to the ground-state. On the basis of this picture 
one may introduce a set of assumptions concerning relative orders of magnitude of certain 
matrix elements, whose validity is obvious in the limit of vanishing two-particle interaction, 



n 



ool 



» KjJ » K/J » (4.7) 

KiJ « Kol, ( 4 -8) 
\n s hI , x \ « KJ if h^I[. (4.9) 

We shall consider diagonal elements to be of zero order, elements connecting states l v to 
I u+P to be of pth order. 

Considering assumption (|4.8| ) first, it asserts that for I belonging to the first few levels of 
the hierarchy, if N, the number of particles is not too small, in lowest approximation matrix 
elements diagonal in / are equal to their value for / = 0. It is easiest to see this for the 
density itself, since the wave functions of the excited states differ from those of the ground 
state by at most a few particles out of N. That this property of matrix elements diagonal 
in / follows for quantities other than the density itself is a consequence of their relation 
to the density, as will be seen from further study below. We shall consider all diagonal 
matrix elements to be zero order quantities. A further assumption, in terms of this scale, 
is that matrix elements in which / and I' belong to adjacent levels in the hierarchy are, 
on the average, of order (1/y/N) compared to zero order quantities. For the sorting of our 
equations, we also need the assumption that matrix elements in which /, I' differ by two 
levels or refer to two different states of the same level are second order quantities, i. e., of 
the order of the product of first order quantities. Of course, it has to be verified a posteriori 
that the solutions found are in accord with these statements. 

Our aim is to apply these assumptions to choose those matrix elements of Eq. ( f4.4|) that 
characterize the state and the states I\. To carry out this program, we must look more 
closely into the structure of the effective interaction \ s . First we rewrite the trace of the 
Hamiltonian in the form 

H = T s + V + W c + H xc , (4.10) 



which defines H xc . It follows that 



\ f dxdx'^n^x)--^— n s ri (x'), (4.11) 



^ = ^ + r + ^ (412) 

= v(x)6 IP + v c IP {*) + i>?£(x), (4.13) 
dx!-. — - — Xr,(x'). (4.14) 



%,(xj 



x — X' 
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The main reason for exhibiting these formulae is to recognize, as we shall see in more detail 
below, that the off-diagonal elements of h are at least linear in the corresponding off-diagonal 
elements of n s . This is obvious from Eq. ( (4.14| ) for the Coulomb contribution and will be 
argued more closely later for v xc . Thus we may safely assume that that the matrix elements 
of h are the same order of magnitude as the corresponding matrix elements of n s . 

Turning then to our major task, which is to study the matrix elements of Eq. (|4.4|) , we 
consider first the ground or 00 element. Neglecting terms of second order and higher, we 
find 

K (xx')^ (x') - /4(xK (xx')] = 0. (4.15) 

It is consistent with our approximations to identify 71q (in leading approximation only) with 
the ground state density of KS theory and /?-oo( x ) with the KS single-particle Hamiltonian. 
Equation (|4.15|) is thus the KS equation in density matrix form and determines a complete 
set of quasiparticle orbitals v ? a(x), where a = h will refer to the quasi-orbitals occupied in 
the ground-state and a = p to those unoccupied. 

Consider next the first-order matrix element 01. Retaining only first-order contributions 
(leading corrections are third order), we may write 

wmS 1 (xx / ) = [n s 00 {xx')h s 01 {x.') + < 1 (xx>* 1 (x') 

-^W^(xx') - ^i(x)« s u (xx')]. (4.16) 

As a first step in the evaluation of this equation, we may, according to Eq. (|4.8|) , set the n\\ 
matrix elements equal to the rioo ones. We also drop the subscripts 00 understanding these 
according to the previous identification to be the standard KS quantities. If we can exhibit 
h s Ql as an (approximate) linear functional of n s Q1 , Eq. (|4.16|) will have the form of a linear 
eigenvalue problem. First we have (the matrix elements in question are local functions of x) 

^ 1 (x) = <(x) + < 1 c (x), (4.17) 

<(x) = / dx'|— L^(x'). (4.18) 

We see that v c is, by definition, already of the desired form. 

We turn then to v xc . Our approach to this quantity is to revert to the study of H xc , 
defined in Eq. ( 4.10 ), which we consider, in line with assumptions previously made, a func- 



tional of noo ~ n, of Uq 1} and of nf , the latter two considered as small quantities. (It is also 
a functional of the other off-diagonal elements, n^, and nl, , where 1' refers to any of the 
other states at level one of the hierarchy of states. It is simply that this dependence does 
not enter into the current discussion.) We then expand H xc as a functional Taylor series in 
these quantities, 

H- = fl-|„ + / (i x^| n; (x) + j^JE-Ln^) 



f 1 5^H XC f 5^H XC 

,1 5 2 H XC 



oi (x' 



^2 ^ 1 (x)^ 1 (x l ^ (x) ^ (x/) + -- (419) 
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Strictly, the quantity H xc \ and its functional derivatives still depend on tin as well as 
n 0Q . It suffices to ignore the difference of the two quantities in the present discussion, but 
we shall have to remember and include the difference in the arguments of Sec. V. We note 
further that only the first and fourth of the terms shown explicitly in this equation are non- 
vanishing. Recall that H xc is a trace and therefore invariant under a unitary transformation 
in the space of states /. Its dependence on the matrix n must also be in the form of traces 
over these indices. As we can see on the example of the Coulomb interaction, this dependence 
is more general than traces of products of n at the same point, but in any event it follows 
that for every factor of n* at some spatial point, there must be a factor of n^, at a generally 
different point. The simplification described above follows. We thus compute to first order 



<(x) 



/ dx faW^M lo " ol(x) 

= J dx.'f WtW (\x- x.'\,n)n s 01 (x.') 

« J dx'/dx-xVKtx'). (4.20) 

In passing from the second to the third line of this equation, i. e., in ignoring the state- 
dependence of /, we are making an approximation equivalent to the adiabatic approximation 
widely used in TDDFT. With the definition (the dependence on n being understood) 

/ e "(|x - x'|) = -J— + /(|x - x'|), (4.21) 
x — X 



Eq. Q4.16|) may be rewritten as (using /i s (xx") = <5(x — x')/i s (x) for convenience) 

um' 01 (xx!) = J dxV(xx')/ e// (|x' - x"IKi(x") + ™oi( xx '>Vx) 

-^(xx")<(x'V) - ^(xx')/ e// (|x - x"|)n;(x")]. (4.22) 



The final task with respect to this equation is to convert it into a standard RPA form. 
Toward this end we reexpress the matrices n s and rigi in terms of the KS single-particle 
functions, y? a ( x ), satisfying the KS equation 

J dxV(xx> a (x') = e a y? a (x). (4.23) 
First of all we have the familiar equation 

n s (xx') = E^(x)^(x'). (4.24) 

h 

Next we must evaluate the sum 

<(xx') = £^( x °V}( x, i)- (4.25) 
j 

Here we must introduce assumptions concerning which values of J contribute to the required 
order. In the space of the eigenstates of the fully interacting system, we are concerned with 
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the ground state and with states that are quasiparticle-hole (qph) excitations of this state. 
When we remove one particle from (create a hole h in) such a state, we expect to encounter 
states that can be characterized as either Oh or 1/t, and these are the values of J that we 



assign in the sum (|4.25|) . If we consistently use the approximations <foh(0) ~ fihO-) ~ Vh, 



the weak-coupling value of Eq. ( [4.25[ ) becomes 

'oi(xx') = EKW4(x'l) + <M*0K(x')]. (4.26) 



Tin 



The final form for this quantity is achieved by expanding the first-order amplitudes in 
terms of KS modes, 

<p 0h (l) = 'Z<p P X p h, ( 4 -27) 
p 

Mo)-E* ( 4 - 28 ) 

p 

The restriction of the sums on the right-hand sides of these equations is also consistent 
with the weak-coupling picture painted above. Strictly the amplitudes X, Y should carry 
superscripts 1, identifying the eigenstate to which they refer, but we shall suppress these 
except when required for clarity, as in Sec. V. Finally then, 

<(xx') = Y& h {*)v;{*!)x; h + ^;(x)^(x')y;j. (4.29) 

p,h 



Introducing Eqs. ( f4.24| ) and ( |4.29j ) into Eq. ( |4.22[ ), we can project out equations for X* h 
and Y* h . We quote the complex conjugate of these equations: 

(e/i - e p + LOi)X ph = (f eff ) p h'hp'X p > h i + (f eff )p P 'hh'Y p 'h', (4.30) 
— t p — LUi)Y ph = (f eff )hp'ph'Y p 'h> + {f eff )hh'pp'Xpi h > 1 (4-31) 

(f eff ) abcd = J rfxrfxV:(x)^(x')/ e// (|x - x'|)^ c (x)^(x'). (4.32) 

The equations found are of the same form as those of the random phase approximation 
(RPA) and agree in detail with the eigenvalue equation that has been derived from the 
density-matrix version of CKS theory. Solutions are to be normalized in the usual way, 
according to the conditions (Appendix B), 

J2(\X P h\ 2 - \Yp h \ 2 ) = 1. (4.33) 

As is well known, two different non-degenerate solutions of the RPA equations are orthogonal 
with the same metric as in ( |4.33| ). 

It is important to emphasize what has been accomplished by the calculations of this 
section. With the help of Eq. ( |4 . 2 9| ) , for instance, we can calculate the off-diagonal matrix 



elements of the density operator between the ground state and the first level of excited states. 
This can be applied, for example to the calculation of the corresponding matrix elements of 
the electric dipole moment. We have yet to establish, however, that the eigenvalues uji can 
be identified as the excitation energies of the system. We turn to this task now. 
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V. EXCITATIONS AS ENERGY DIFFERENCES 



In principle the energy differences can be calculated from the expression 



= J2(I\H\I)-2{0\H\0) 
1=0,1 

= Ei — E , 



(5.1) 



where Ej is the energy of state I. This difference will be evaluated with the aid of Eqs. ( |4. 10[ ) , 
( |4.11| ), and the simplified version of (|4.19| ). The result that we shall establish is 

Ei — Eq — ^{(e p — €h)(\X ph \ 2 — |l^h| 2 ) + X* h [f ph i hp iX p i h i + fppfhh'Yp'h'] 

p,h 

+Y p *h[fhp'ph'Ypih> + fhh'pp'Xpih']}- (5.2) 
But the right hand side of this equation is easily seen from Eqs. ( |4.30| ) and ( |4.31| ) to equal 



toi , provided that we make use of Eq. ( 4.33 



It is simplest to evaluate the difference ( |5.1| ) first for the interaction terms. Consider, for 
instance, the Coulomb difference, 



AW C 



J dxdx'--, 



2 x 



[n^fxln^fx' 



ng (x)nnn(x') 



2< 1 (x)^ (x')] 



/ ^ x/ | x ^ x ,| {Ki(x) ~ ^oo(x)Ko( x ') +<(x)^ (x / )}, 
J &K 1 (x)-^ (x)]t> c (x) + J rfxrfx' ^ - ng 1 (x)^ (x / )]. 



(5.3) 



To obtain the value exhibited in the first line, we have made the quasi-boson approximation 
n s 12 = y^^oi) which is an expression of the closure approximation referred to in Sec. II and 
discussed in more detail below. The further simplification is made possible by the fact that 
the difference n\ 1 — n s m (see below) is quadratic in the RPA amplitudes. The corresponding 
difference involving the exchange-correlation energy can be written 



rfxK^x) - r^ (x)K c (x) + / dxdx7(|x - x'|)n^ 1 (x)^ (x / )], 



AW* 



Next we see that the second terms of Eqs. (|5.3|) and ( |5.4|) combine to give 



(5.4) 



rfxc/x7 e// (|x-x'|)n^( X ; 



n io l x 



Xphifph'hp'Xp'h' + fpp'hh'Yp'h' 
J T~Yph[fhp'ph'Yp'h' + fhh'pp'Xpi, 



(5.5) 



which has been evaluated with the help of Eq. ( |4.29| ). This is already seen to be the 
interaction terms of Eq. ( |5.2|) . 

The remaining terms of Eqs. ( |5.3j ) and ( |5.4| ), as well as the contributions arising from 
the kinetic energy and the external potential depend on the value of 



%l(XJ 



oo 



x ) = £b}( xl VXxl) - ^(xO)^(xO)]. 



(5.6) 
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To enumerate the states J that contribute to this difference it is useful to picture the state 
1 as an elementary boson excitation, as is done in the standard approach to the RPA. The 
relations that follow from this assumption will lead, as we shall see, to a quantitative form 
of closure approximation that is essential to the calculation. By the notation 1 x 1, we 
shall mean a double boson excitation with the same boson, whereas by 1 x 1' we shall mean 
a double excitation with different bosons. Thus for the amplitudes y?j(l), we consider the 
values J = Oh, lh, 1 x lh, 1 x l'h. The contributions from the latter two choices are evaluated 
in boson (closure) approximation as 

(1) = v^iCO), (5.7) 
¥>ixi'h(l) = ¥>i'(0)- ( 5 - 8 ) 



For the amplitude <fj(0), the required values are J = Oh, lh, l'h. For the difference ( f>.6\) , 
we thus find, suppressing coordinate indices, 

n s n ~ noo = £[<^(1 Wl) + vlfc(OWO) + ^(1)^(1) - <pM<Poh(0)]- (5.9) 



The total contribution of the first two terms of Eq. (|5.9|) to the energy difference un- 
der study, obtained by substituting Eqs. ( [4.27]) and (|4.28| ) and applying the result to the 



sum of single-particle operators that add up to the KS Hamiltonian h s , is found to be 
J2 p ,h[ e p{\X p h\ 2 + |^p/i| 2 )], one of the single-particle terms in Eq.(|5.2|). The evaluation of the 
remaining terms of Eq. (|5.9|) is carried through by studying the normalization conditions, 
Eq. ( |3.11D . We calculate 



£W)f 



i 



l</Wo)| 2 + l</Wi)l 2 + £ l</Wi')l 2 , (5.io) 
£W)I 



2 



1 



Wih{l)? + \Vih(0)\ 2 + \<p lh (l x 1)| 2 + £ \<p lh (l x l')| 2 

|^(1)| 2 + |^(0)| 2 + 2|^(1)| 2 + £ l^ohCl')! 2 . (5-11) 



where the last evaluation has made use of the boson approximation expressed by Eqs. 
and (|5.8|) . These equations are satisfied by introducing a renormalization of the KS orbitals 

m (x0) = ^(x)[l -W\X p h\ 2 - ^£ £ l^il 2 ], (5-12) 

^(xi) = ^(x)[i - Yl l^l 2 - \ £ I VI 2 - \ £ £ l^il 2 ]- (5.i3) 

p z p z p V ^x 

Combining these results and applying them to the last two terms of Eq. ( |5.9j ), suit- 
ably multiplied by the sum of terms that comprise h s leads to the final contribution 
Y,ph[-£h{\X ph \ 2 + | VI 2 )] to the theorem stated in Eq. Q. 
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VI. CONCLUDING REMARKS 



In this paper, we have developed an alternative formalism for the study of excited states 
within a framework that generalizes the basic ideas of KS theory. The main novelty in our 
approach compared to other methods is that the latter work with a single density, be it the 
average in the ground state, in an excited state, an ensemble average, or the average in a 
suitably chosen time-dependent state. On the other hand, we arrive at a formalism involving 
an entire array of matrix elements of the density operator taken among a pre-selected set 
of states. The application of the variational principle for the trace of the Hamiltonian then 
leads to a generalized KS scheme in terms of orbitals that depend not only on the coordinate 
x, but also on a label I of the included states. We have examined the consequences of 
this formalism for the weak-coupling limit. We did this by framing a set of assumptions, 
including a closure approximation, in order to identify the most important amplitudes and 
their equations that characterize the ground state and a simple class of excited states that 
are composed of lp-lh excitations of the ground state. 

In this way, we regained first the ground-state KS theory and second derived an eigen- 
value equation of RPA form. By approximating a state-dependent (frequency-dependent) 
effective interaction by a state-independent (frequency independent) effective interaction, the 
eigenvalue equation became identical to one that was first derived from the density-matrix 
version of CKS theory [JT2 . 



In our formalism, it is not immediately obvious that the eigenvalues, which originally 
entered the theory as Lagrange multipliers in a variational principle, can be identified with 
excitation energies of the physical system. We prove that this is correct identification, in 
agreement with previous results. The theory also provides the means for the computation of 
electromagnetic transition rates from the excited states to the ground state. As formulated, 
the theory described in this paper can be extended to improve the approximations made 
for lqp — lqh states, as well as to study more complicated states, for example, of 2qp — 2qh 
character. 
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APPENDIX A: RELATION OF WEAK-COUPLING LIMIT TO 
TIME-DEPENDENT DENSITY FUNCTIONAL THEORY 



In this section, we shall connect the linearized RPA equations ( |4.30|) and ([4.31|) with 



a corresponding linearized approximation to TDDFT. We start with TDDFT in density- 
matrix form 

^ = [(r + ^)),A (Al) 

p s (xt,x't) = ^^(xt)^(x't), (A2) 
h 
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v'(xt) = sAjj(V(t) + W(t) + T(t) - T s (t)). (A3) 

Here <p(xt) are the N instantaneous eigenfunctions of r + v s (t) of lowest energy, defining a 
time-dependent Slater determinant whose kinetic energy is T s (t), and V(t), for example, is 
the expectation value of V in the time-dependent wave-function 

We are interested in the physical situation where the time-dependence of the state vector 
arises not from an explicitly time- dependent external field but from the fact that initially the 
state vector is a superposition of the ground state (predominately) and a small amplitude 
for one of the excited states. We thus assume that 

p s (xt,x't) = p°(x,x') + [p 1 (x,x / )exp(-?cjt) + ex.], (A4) 

^(x,xO = £[AVy? p (xK(x') + y p ^(xv;(x')]. (A5) 

ph 

In (|A4j) and below the superscript identifies quantities associated with the KS ground-state 
theory. 

What follows now is close to a standard derivation of the RPA. We insert (|A4j) and (|A5|) 
into (|A1|) and, considering the amplitudes X and Y as first order quantities, we expand to 
first order. For this purpose, we need the expansion, 

v s (-xt) = u°(x) + J dx'/(x,x')n 1 (x') J (A6) 

= S?g- (A7) 

n 1 (x)=p 1 (x,x). (A8) 

In Eqs. ( |A6p and (|X7|), we have already made the adiabatic approximation by ignoring 
the time dependence of /. As a consequence, the quantity called / in this appendix can 
be identified with the quantity f e ^ of the text. From the zero order term, we regain the 
KS theory for the ground state. From the first order terms proportional to exp(— iut), for 
example, we find 

r fiv 

cV(x,x') = [(r + ApW') + J dx"[j^? ^-yP'K^W). (A9) 

Taking, in turn, the qph and qhp matrix elements of ( |A9| ) , we find the familiar equations 
[e/i — e P + w\X p h = fph'hp'Xpih' + fpp'hh'Yp'h', (A10) 

[ e /i ~ e p ~ u ]Yph = fhp'ph'Yp'h' + fhh'pp'Xpih' ■ (All) 



APPENDIX B: RPA NORMALIZATION CONDITION 

We define mode operators, a a , for the field V>(x) by expanding in terms of the KS modes, 

V^( X ) = H a a<Pa(x), (Bl) 
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a = {h,p}. From the commutation relations for quasiparticle-quasihole pairs, 



[a[a p , a p ,a h >] — 5 hh >5 pp ' — S hh >a p/ a p — 8 pp >a h 'a h , (B2) 

we obtain an approximate sum rule by taking the expectation value in the state |0), intro- 
ducing a complete set of intermediate states and retaining only the first term on the 
right hand side (on the justified assumption that, for instance, (0| a^a^|0) is, on the average 
small compared to unity). With the definitions 

4=<0|4a p |z), (B3) 

V; h = (0\ala h \t), (B4) 

we have 

E [£h#h> - ih'Vpl) = V<W ■ (B5) 

i 

We would like to identify the quantities £ and 77 with the quantities X and Y, where the 
latter satisfy Eqs. ( |4.30| ) and (|4.31|) . Equation ( |B5| ) would then constitute the completeness 
relation for the solutions of these equations, and as is well-known, a completeness relation and 
orthogonality of solutions with the corresponding metric implies the normalization condition 
Eq. ( ^4.33|) . Toward this end, we consider two different evaluations of (O|^'''(x)'0(x)|z) = 
n,o(x). On the one hand we have in an approximate evaluation based on the physical 
picture, 

^o(x) = E^a( x V&( x )(0|ala&K) 



ah 



EK( x )^( x K°K a ^)] + ^( x V P ( x )(°l a l a pK>- (B6) 

ph 

On the other hand, from the generalized KS mapping n«o — > n s i0 and Eq. ( |4.29| , we have 

n, (x) = ^[^(x)^(x)y;, + ^(x)^(x)Xy. (B7) 

ph 

The identifications £ = X and rj = Y are consistent with these equations. 
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